
// $Id: HJM_Main.cpp 2 2012-05-22 17:55:16Z gerry@gaboury.biz $

#include "Matrix.hpp"
#include "Util.hpp"
#include "HJM.hpp"
#include "HJM_Observer.hpp"
#include "Config.hpp"
#include <stdexcept>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <unistd.h>

// Boost
#include <boost/thread.hpp>
#include <boost/bind.hpp>
#include <boost/shared_ptr.hpp>

using namespace std;

int
main ( int argc, char* argv[], char* env[] )
{
  // defaults
  int nPath = 1;
  int nStep = 1000;
  double Tf = 10.0;
  ulong seed = 0;
  int nThreads = 0;
  int c;

  while ((c = getopt(argc, argv, "f:n:p:s:t:")) != EOF)
    switch (c)
      { 
      case 'f':
        Tf = atof ( optarg );
        break;
      case 'n':
        nStep = atoi ( optarg );
        break;
      case 'p':
        nPath = atoi ( optarg );
        break;
      case 's':
        seed = atol ( optarg );
        break;
      case 't':
        nThreads = atoi ( optarg );
        break;
      case '?':
        cerr << "invalid option" << endl;
        exit(1);
      }

  if ( argc - optind < 2 )
    {
      cerr << "Usage: " << argv[0] << " <initial rates - F0.txt> <input file>" << endl;
      return 1;
    }

  try
    {
      //
      // If seed not initialized, use time.
      //
      if ( ! seed )
	seed = static_cast<ulong>(std::time(0));

      //
      // setup threading
      //
      if ( ! nThreads )
	nThreads = boost::thread::hardware_concurrency();
      cout << nThreads << " threads (cores)" << endl;

      vector<boost::shared_ptr<boost::thread> > m_thread;

      // what do you want to price?
      Config* CFG = Config::instance();
      CFG->read(argv[optind+1]);

      //
      // read input data
      //
      ifstream HJM_data("data/HJM.dat",ios::binary);
      if ( HJM_data.fail() )
	{
	  cerr << "open \"data/HJM.dat\" failed" << endl;
	  exit(1);
	}
      int nX;
      HJM_data.read((char*)&nX,sizeof(int));
      cout << "Reading data for " << nX << " maturities" << endl;
      double Vtmp[nX];
      // read maturities
      HJM_data.read((char*)Vtmp,nX*sizeof(double));
      Vector<double> X(nX,Vtmp);
      // read drift
      HJM_data.read((char*)Vtmp,nX*sizeof(double));
      Vector<double> M(nX,Vtmp);
      // read eigenvectors
      Matrix<double> EigenVec(3,nX);
      for ( int i=0; i < 3; ++i )
	{
	  HJM_data.read((char*)Vtmp,nX*sizeof(double));
	  EigenVec[i] = Vector<double>(nX,Vtmp);
	}
      HJM_data.close();

      // initial rates
      ifstream F0_data(argv[optind]);
      if ( F0_data.fail() )
       	{
       	  cerr << "open " << argv[optind] << " failed" << endl;
       	  exit(1);
       	}
      Vector<double> F0(nX);
      F0_data >> F0;
      F0_data.close();

      //
      //  check that everything was read correctly
      //
      ofstream ofs_vol("vol1.txt");
      if ( ofs_vol.fail() )
	{
	  cerr << "open \"vol1.txt\" failed" << endl;
	  exit(1);
	}
      ofs_vol << "Maturity";
      for ( int j = 0; j < nX; ++j )
	ofs_vol << " " << X[j];
      ofs_vol << endl;
      for ( int i=0; i < 3; ++i )
	{
	  ofs_vol << "Vol_" << i;
	  for ( int j = 0; j < nX; ++j )
	    ofs_vol << " " << EigenVec[i][j];
	  ofs_vol << endl;
	}
      ofs_vol << "M";
      for ( int j = 0; j < nX; ++j )
	ofs_vol << " " << M[j];
      ofs_vol << endl;
      ofs_vol.close();

      // 
      // HJM simulation itself
      //

      // Initialize the HJM simulation
      int nPathsPerThread = nPath / nThreads;
      if ( nPathsPerThread <= 0 )
	nPathsPerThread = 1;
      cout << "\n*** " << nPathsPerThread << " simulations per thread ***\n" << endl;

      HJM* iHJM = new HJM(X,M,EigenVec,F0,Tf,nStep,nPathsPerThread,nThreads);

      // Initialize the HJM observers and attach to the HJM object
      vector<boost::shared_ptr<HJM_Observer> > Observers;
      while ( CFG->next() )
	Observers.push_back(boost::shared_ptr<HJM_Observer>(new HJM_Observer(iHJM,CFG->element())));

      Matrix<double> RateMean;

      for ( int ith = 0; ith < nThreads; ++ith )
	{
	  m_thread.push_back(boost::shared_ptr<boost::thread>
			     (new boost::thread(boost::bind(&HJM::simulate,iHJM,ith,seed))));

	  //Matrix<double> Rates = iHJM->result();
	  //RateMean += Rates;
	}
      // Join threads
      for ( int ith = 0; ith < nThreads; ++ith )
	m_thread[ith] ->join();
      

      //vector<double> TS = iHJM->TimeSeries();

      //RateMean /= nPath;

      //
      // output pricing/option results
      //
      for ( int i=0; i < Observers.size(); ++i )
	{
	  vector<double> args = Observers[i]->args();
	  cout << Observers[i]->name() << "(";
	  cout << args[0];
	  for ( int j=1; j < args.size(); ++j )
	    cout << "," << args[j];
	  cout << ") = " << Observers[i]->value() << endl;
	}

      // ofstream ofs_out("HJM_out.txt");
      // if ( ofs_out.fail() )
      // 	{
      // 	  cerr << "open \"sim.txt\" failed" << endl;
      // 	  exit(1);
      // 	}
      // ofs_out << "Term";
      // for ( int j=0; j < nX; ++j )
      // 	ofs_out << " " << X[j];
      // ofs_out << endl;
      // for ( int i=0; i < nStep+1; ++i )
      // 	{
      // 	  ofs_out << TS[i];
      // 	  for ( int j=0; j < nX; ++j )
      // 	    ofs_out << " " << RateMean[i][j];
      // 	  ofs_out << endl;
      // 	}
      // ofs_out.close();

      delete iHJM;
  }
catch ( std::exception& ex )
  {
    cerr << ex.what() << endl;
  }
  return 0;
}
